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I. INTRODUCTION 

Progress toward reliable and efficient finite element 
capabilities for thin shell instability analysis can at best 
be described as uneven. A major reason for the slow rate of 
accomplishment is that the problem encompasses numerous dis- 
tinct facets, each of which comprises in itself an exceptional- 
ly complex development challenge. The three major facets of 
finite element instability analyses are: 

a) representation of shell geometry 

b) representation of shell deformational behavior, and 

c) algorithmic tools for solution of the large-order 
systems of nonlinear equations characterizing the 
various phases of shell instability. 

These and other aspects of the problem have been the subject of 
the grant research of which the present report forms a part; 
the specific goal of the research being a working computational 
capability for the elastic instability analysis stiffened thin 
shells of rather general configuration. 

The present report is concerned with the formulation of 
the stiffness relationships for a doubly-curved triangular thin 
shell element. Of necessity, it deals with both the representa- 
tion of shell geometry as embodied in the element and with the 
representation of deformational behavior, i.e., with items (a) 
and (b) , above. It is pertinent to note that algorithmic tools 
for the present research (item (c)) have been described in two 
previous NASA Contractor Reports (Refs. 1 and 2); a number of 
published papers (Refs. 3-5) and further publications on this 
aspect will appear in the future. 

Formulation of total finite element representation for 
geometrically nonlinear analysis is felt to be of such scope 
that two reports are needed for adequate coverage of the mate- 
rial v The present report deals only with formulative aspects 
for linear analysis for a triangular shell element. Nonlinear 
and instability aspects are treated in a companion report (Ref. 6). 
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Some preliminary comments are in order with respect to 
the developmental history of finite element representations 
for thin shell analysis, since this background is the predom- 
inant influence, upon the direction taken here. One of the 
first surveys of finite element representations for thin shell 
analysis was published in 1969 (Ref. 7). A significant number 
of other surveys (Refs. 8-10) have since appeared. These dem- 
onstrate that nearly all researchers have preferred to follow 
the well-established path of minimum potential energy (dis- 
placement-based) formulations and that extreme difficulty is 
encountered in attempting to satisfy the necessary compati- 
bility conditions. With this is mind the present work de- 
parts somewhat from the conventional approach to element for- 
mulation and adopts a "generalized" potential energy scheme 
(Refs. 11 and 12). The resulting formulation is relatively 
simple and in the performance of numerical analyses is shown 
to be of acceptable accuracy. 

The report is organized as follows. A review is first 
given of alternative formulations for thin shell finite ele- 
ments. Because of the extensive information available in 
Refs. 7-11, this review focuses only on the more recent and 
sophisticated developments. The review establishes further 
motivation for the approach adopted here. The geometry of 
the subject shell element is then defined and succeeding 
sections examine the underlying shell strain-displacement 
equations, the variational principle employed, the chosen 
displacement functions, and special aspects of the solution 
process. Rather extensive numerical solutions and compari- 
sons with alternative formulations are presented in the con- 
cluding section. 

II. ALTERNATIVE FORMULATIONS 

A formidable literature dealing with the establishment of 
finite elements for thin shell structural analysis has emerged 
since 1967. The aforementioned review papers have served to 
classify the respective formulations and to identify their 
relative merit. The present section is devoted to an examina- 



3 


tion of the more sophisticated and accurate of the formulations 
available in the literature at the present time, as identified 
in these reviews or identified as such by application of criteria 
given in these reviews. Three formulations, each of triangular 
shape, and each based on assumed displacement fields and the 
principle of stationary potential energy, are described. 


The first of the curved thin shell finite elements discussed 
in this section is due to Cowper, Lindberg, and Olson (Refs. 13, 
14) . Originally formulated in terms of shallow shell theory 
(Ref. 13) , this element was extended (Ref. 14) to cover nonshallow 
applications. A 21 -term quintic polynomial was chosen for the 
normal displacement field, w, and a 10-term cubic field for each 
of the in-plane displacements, u and v. By imposing cubic normal 
rotations along each edge, the derivation satisfies interelement 
continuity. The centroidal values of u and v can be solved for 
and eliminated from the stiffness matrix. The resulting stiff- 
ness matrix then has ■ 36 degrees of freedom, 12 at each of the 


corner nodes, consisting of u, u , u , v, v , v , w, w , w 


XX * 


w and w 
xy yy 


The choice of displacement fields was motivated by the view 
that the consistency of the order of the error term in the strain 
energy expression for a general displacement configuration is 
more important than achieving a zero error for one particular 
displacement mode. Because the second differential of w and only 
the first differential of u and v are present in the strain 
energy expression, the expansion for w is chosen to be one order 
higher than that for u and v. No explicit consideration has 
been made in the derivation for rigid body modes of displacement. 
Note that although a quintic field was chosen initially for w, 
because of the constraints imposed on it, the function can only 
be thought of as a complete quartic with a few extra quintic 
terms . 


The second element discussed in this section is the so- 
called •SHEBA' element formulated by Argyris and Scharpf 
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(Ref. 15). Complete 21-term quintic polynomials are utilized 
for all three displacement components. Interelement continu- 
ity is assured by the use of midside nodes having all three 
normal derivatives as degrees of freedom. The corner nodes 
have 18 degrees of freedom each: 3 translational displacement, 
all six first derivatives and all 9 second derivatives. The 
condition of zero strain under rigid body motion is satisfied 
exactly by using the same interpolating procedure for the sur- 
face location as for the displacements. 

Published results (Ref. 16) demonstrate that 1 SHEBA 1 
produces extremely accurate solutions, but in a recent com- 
parison study (Ref. 17) the plate bending form of this element 
(known as TUBA) proves to be no more efficient than most other 
plate bending elements available, despite being more accurate 
for a given grid refinement. The reason for this is that the 
bandwidth of the algebraic equations of the problem, when ex- 
pressed as a percentage of the total number of structural 
freedoms, is very high for an idealization using TUBA (or 
SHEBA) elements, in comparison to simpler elements, and the 
effort required to solve a system of equations is proportion- 
al to the square of the bandwidth. One must also take into 
account the effort required to construct the individual ele- 
ment stiffness coefficients. 

The third and final triangular thin shell finite element 
to be commented upon here has been published by Dupuis and 
Goel (Ref. 18). This formulation satisfies all necessary 
requirements for convergence of the potential energy. Ra- 
tional functions are used for the displacement functions as 
well as to define the position of the shell above a datum 
plane. The use of rational functions other than simple poly- 
nomials is described in the text by Zienkiewicz (Ref. 20) 
where they are called singular functions. In this element 
the rigid body motion condition is satisfied in the same way 
as for the SHEBA element above. Two levels of sophistication 
of the element are available, the first having the required 
continuity of first derivatives, whereas the second has 
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continuity of second derivatives, which are of course exces- 
sive for the imposition of the variational principle. Dupuis 
and Goel conclude that the second level i$ the most efficient. 

In a later paper Dupuis (Ref. 19) has used the same ap- 
proach to develop an almost identical element. The latter 
derivation utilizes cubic rational functions which are some- 
what simpler than the previous ones. The element has nine 
degrees of freedom at each of the corner nodes. Comparison 
results are given for a cylindrical shell test problem. 

Each of the above elements is complex from the stand- 
point of formulative effort. This is a disadvantage in ex- 
tension to nonlinear problems. The times required for gener- 
ating the nonlinear stiffness terms are higher than for the 
linear terms (assuming the same displacement fields are used 
in calculating both the linear and nonlinear terms) and may 
become prohibitive. Thus, we seek in the following a curved 
triangular thin shell element formulation which preserves 
simplicity without major reduction in accuracy. Before this 
development is presented, however, we review briefly in the 
next section the mode of description of the element geometry. 

I II. GEOME TR IC RE PRESE NTATION OF ELEMENT 

Investigators interested in shallow shells have often 
made the assumption that the shape of each element in the 
shell idealization is defined by the projection of the ele- 
ment on the base plane. The curvatures are either input di- 
rectly or are calculated from information about the height 
of the shell above the base plane. This approach has the 
disadvantage that small gaps exist between elements for a 
doubly curved shell. Also, it cannot be extended to cover 
deep shells because of the possibility of elements being per- 
pendicular to the base plane. 

Lien (Ref. 21) has used Coon's (Ref. 22) idea of surface 
patches to define the geometry of the shell. By selecting 
the quadrilateral patches to represent elements, and using 
the same bicubic functions for the displacements and geome- 
try, a powerful tool is developed. 
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In the derivation discussed in this report, the curva- 
tures are usually specified directly, because for most 
shells encountered in practice (translational and rotation* 
al) , these quantities are easily evaluated. If, for a par- 
ticular structure, they are not, the surface patch approach 
of Ref. 21 can be used as an option. 

The basic geometry of the triangular shell element is 
shown in Fig. 1. Points are located on the shell surface 
in the orthogonal curvilinear system a-0. The user provides 
nodal values for the curvatures 1/R^, I/R 2 an< * as well 

as the shell thickness t. These variables are interpolated 
linearly over the area of the element. Because numerical 
integration is used in the stiffness matrix calculation, any 
other interpolating function could be used. 

IV. SHELL THEORY 

As in any particular aspect of structural mechanics, 
the governing differential equations of thin shell analysis 
are derived by combination. of the following component rela- 
tionships 

i) Strain-displacement 

ii) Equilibrium 

iii) Stress-strain 

The proper forms of relationships i) and ii) have been 
the subject of disagreement among shell theorists for some 
time and still other disagreements arise on account of ap- 
proximations made after combination of the above to form 
the governing differential equations. Nevertheless, due to 
the work of Koiter (Ref. 24) and others, an understanding 
has been reached of the significance of these differences. 
Koiter has shown that many investigators have differed in 
their retention of terms of the order of magnitude of e/R, 
where e represents a membrane strain component and R is the 
corresponding radius of curvature of the shell. Despite 
this, their seemingly contradictory results are equivalent 
in the context of linear thin shell theory whose basic 
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approximations give rise to an error in the strain energy 
of the same order of magnitude. This error is caused in 
the main by the neglect of the strain energy corresponding, 
to the transverse shear and direct stresses. Some writers 
have retained terms smaller than e/R while adhering to the 
basic approximations of thin shell theory (the Love-Kirch- 
hoff postulates) . This inconsistent retention of smaller 
terms may well produce improved accuracy for particular prob 
lems, but for the general case their retention cannot be 
justified. 

The present development deals with the appropriate form 
of the strain energy, rather than the form of the governing 
differential equations. There is general agreement on the 
expression for strain energy for linear thin shell analysis, 
which can be written as follows. 


U 


1 

1 

1 

1 


— + e 2 ) 2 - 2(1 - v)( ei e 2 - \ * 2 )] 

Et ~ n [(*. + K .) 2 - 2(1 - v)(k.k 2 - t 2 )] 

12(1 - v*) 1 z 1 L 


( 1 ) 

where the in-plane strains e and <f> , and the curvatures k 
and x are defined in Eq. (2) , below. 


To transform Eq. (1) into the form needed for a finite 
element formulation we must introduce the strain displace- 
ment equations. As noted above, a variety of different 
forms of these equations have been proposed for linear thin 
shell analysis. The strain-displacement equations given by 
Koiter (Ref. 24) are used in the work reported herein. 
Koiter's relationships are described as "consistent”, which 
means that the expressions are neither deficient in terms 
larger than e/R, nor do they contain superfluous terms small 
er than e/R. Thus a consistency with the accuracy of the 
Love-Kirchhof f postulates is achieved. These relationships 
reduce exactly to those presented by Novozhilov (Ref. 25) 
when the coordinate system is chosen to coincide with the 
principal curvature directions (i.e. when 1/T - 0). 
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The relevant relationships are: 


1 3u v 3A _ w 

A Fa XF FF " FT 


1 3v + u 


F FF 


3B 

XF FcT 


w 


♦ 

2 n 
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1 3v 1 3u 

X Fa F FF 


u 3B 

XF Fa 


v 3B 

XF Fa 

v 3B 


1 3v . 1 3u . u 3A + 5 

X Fa FFF XF FF XF Fa 




1 3w . u . v 

X FS iq T 


1 3w 

F FF 


1 3 *1 

x Fc r 

, 3*o 


♦ v ♦ u 

*r 

*2 3A . ft 

XF FF T 


$ 


1 3B 


ft 


2 " F FF” + XF Fa ’ T 


3 * 




2w 

T 


? T « 1 2 . 1 W "1 "1 3A 2 3B ,1 1 . r7 . 

Zt XFT" FFF“"XFFF‘XFFS’ (2) 

where A and B are defined such that Ada and BdB are incre- 
ments of length in the a and B directions, respectively. 


It is essential in the intended application that the 
chosen strain-displacement equations do not yield strains 
under the imposition of rigid body motion. Koiter's theory 
is satisfactory with regard to this point. Other theories 
which meet this condition as well as some which do not, are 
identified in Ref. 26. 


One assumption that has been widely used in finite ele 
ment shell analyses is that of shallowness of the shell. 
Resulting from this is the simplification in the curvature 
terms of the strain displacement relationships, as follows: 
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i a^w 

3a 2 

i a 2 w 

b 7 as 7 

i a 2 w 


(3) 


* AB 3a3 B 

By comparison with the appropriate equations in (2) we can 
see that the simplifications are justified if 

1 3 2 w 1 au 

A 7 9a 7 So 


and 


1 3 2 w 

B 7 3B 7 


>> 


i av- 

RJ 38 


i a 2 w >> L_ iH + 1 3v 
AB 3a36 R 3 " 36 R 2 3a 


The left hand side of each expression is a curvature 
term, while the right hand side is proportional to the in- 
plane strains. The three conditions will be satisfied if 
bending action predominates over the membrane behavior. This 
is the usual behavior of shallow shells, although it is easy 
to imagine cases where the loading and boundary conditions 
are such that this is not so. It is interesting to note 
that for the popular test case problem of a shallow cylin- 
drical shell roof under dead weight loading, the difference 
in the solutions for key displacement parameters obtained by 
deep and shallow shell theories, respectively, is 3%. Thus 
this case does not give an adequate test of the suitability 
of the respective theories. 

Assumptions are often made about the shallowness of a 
particular element to facilitate geometric description. An 
example of this type of assumption would be the one that ap- 
proximates the area of the curved element with the area of 
the flat plate joining the nodes. This type of assumption 
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improves with grid refinement. This is in contrast to the 
assumptions concerning the relative magnitudes of the bend- 
ing and in-plane strains quoted above that do not improve 
with grid refinement. Hence, although an element may be 
considered to be shallow, it may form part of a deep shell, 
and therefore the appropriate general strain-displacement 
relationships should be used. This inclusion dictated the 
choice of Koiter's expressions in the present study. 

V. BEHAVIOR REPRESENTATION 

The same displacement functions have been chosen for u, 
v, and w (see Figure 1). These functions are expressed in 
terms of the triangular coordinate system L^, L 3 (Figure 
2). For the u-displacement we have 



u - |Nj {u} 

(4) 

where 



l_Nj * 

L N 1 N al N B1 N 2 N a2 N 32 N 3 N a3 N B3 N 4J 

(4a) 

(u} T * L U 1 

f 3Ux r 3Ux r 3u. r 3Us „ f 3Us f 3u. 

^3x^ 1 ^3g^ 1 u 2 ^3a^ 2 ^33^ 2 u 3 ^3(1^33^ u 4-* 

(4b) 


Note that the subscript 4 refers to the centroidal node 
and the subscripted a f s and 0's denote the nodal coordinates 
in the curvilinear surface system 

♦ 3 L 2 l 2 ♦ SL^Lj - 

N g i « -a 3 (L 2 L 2 -L 1 L 2 Lj) + a^LjLj-l^LjLj) 

N ol “ b 3 (L 1 L 2' L 1 L 2 L 3 ) " b 2 (L l L 3‘ L l L 2 L 3 ) 

where * (flj-Bi) and b 2 * (a^-dj) 

and likewise for N 0 2 » e *c. by cyclic rotation of 

the subscripts. 


Also 


N 4 - 27 44L3 
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The expressions for v and w are of the same form as u, i.e. 

v “ l N j (4c) 

w - jNj {w} (4d) 

where u Nj in both cases is given by Equation 4a and {v> and 
{w} are similar in form to {u> as given by Equation (4b). 

Each element stiffness matrix is formed in the manner 
described in the next section, using the above. 

The reasons for this choice of displacement field are 
as follows: 

1) The cubic function was felt to be the simplest rep- 
resentation with acceptable solution accuracy for the antici- 
pated applications. In theory a quadratic function may be 
used with the flexural strain energy terms (the second deriva- 
tives) and a simpler element is produced (Ref. 27) , but at 
the price of a slower rate of convergence. 

2) Functions of the same order must be used for bend- 
ing (out -of -plane) and axial (in-plane) displacements in the 
modeling of an eccentrically stiffened shell in order to 
prevent relative displacements between the stiffener and 
shell along the line of contact. Retention of functions of 
the same order for both out -of -plane and in-plane displace- 
ments insures interelement continuity of u, v and w displace- 
ment components, even if there is a discrete change of slope 
between elements. 

3) A majority of shells encountered in practice behave 
in a predominantly membrane mode. Thus it would seem that 
reducing the error in the stretching energy is more profit- 
able than reducing the bending energy error by a correspond- 
ing amount. The stretching energy involves lower order 
derivatives than the bending energy so that use of the same 
order of function in approximation of both energies yields 

a higher-order representation of stretching energy. Con- 
clusions drawn by Fried (Ref. 28), using a rigorous approach 
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to error analysis in shell problems, are in agreement with 
the above philosophy. 

4) Mallett (Ref. 29) has shown, in conjunction with 
geometrically nonlinear problems, that the Euler equations 
of the associated functional imply a higher-order represen- 
tation for membrane action than for bending behavior. This 
effect is in contrast to linear analysis, represented by 
item 3, above, but also emphasizes the desirability of mem- 
brane functions of order no lower than the bending functions. 

Note that the rigid body modes that cause no straining 
have not been included explicitly in the displacement fields. 
A great deal of controversy has raged about the importance 
of such modes for curved shell elements. It has been pointed 
out (Ref . 27) that elements that do not meet the condition of 
zero strain under rigid body motion violate equilibrium. 
Convergence rates have been shown to improve when this re- 
quirement is met (Ref. 30 ) . 

It is possible to estimate quantitatively how well ele- 
ments approximate the rigid body modes by evaluating the six 
lowest eigenvalues of the element stiffness matrix (Ref 
13). In general, if higher order expansions are used for 
inplane as well as bending displacements, six of the lowest 
eigenvalues have values sufficiently small to be taken as 
zero for practical purposes. This has been verified for 
axisyrametric shells (Ref. 32) . An attempt is made to ex- 
plain this fact analytically in Appendix A. 

VI. VARIATIONAL PRINCIPLE AND CONSTRUCTION OF STIFFNESS 

EQUATIONS 

The principle of minimum potential energy is used as 
the underlying basis of the element derivation in this report. 
It is well known that the continuity requirements demanded 
by this variational principle can be troublesome for the 
case of triangular plate bending or shell elements. Be- 
cause of this some researchers (e.g. Ref. 27) have neglected 
these requirements while others (e.g. Ref. 33 ) have found 
that better accuracy can be expected if the requirements 
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are relaxed. In the present development steps have been 
taken to meet the continuity requirements so that conver- 
gence is ensured for all mesh patterns. 

As Eq. (2) discloses, the second derivative of the 
normal displacement is the highest order derivative to ap- 
pear in the strain displacement equations and hence in the 
strain energy expression. Since a variational principle 
requires that the functional be integrable over the com- 
plete domain, continuity of the first derivative of normal 
displacement is needed to enable the integration of the 
strain energy over the whole shell. 

For ease of discussion, the interelement continuity 
of a special case of the shell element - a flat plate - is 
considered. The same considerations will be valid for a 
curved surface. 


The displacements along the edge of an element vary as 
a cubic function which is uniquely defined by the end values 
and slopes. The rotation of the element about the edge fol- 
lows a parabolic distribution. To define this parabola, 
only two end values are available, and thus displacements 
of nodal points not on the edge are needed for the complete 
parabola definition. Thus at a common edge the parabolas 
for adjacent elements will in general be different. This 
difficulty is overcome by forcing the normal rotations at 
the midside to be the same for adjacent elements. The 
Lagrangian multiplier technique is used to impose these 
constraints (Refs. 11, 34). The necessary modifications 
to the principle of potential energy is outlined below. 


The potential energy for the system, II, is 

n - i L AjtK]{A} - L Aj{P} - f pwdA (S) 

* A 

where (A) is the vector of nodal displacements 
{P} is the vector of applied joint forces 
| pwdA is the contribution to the energy from the 
^ distributed loading 

[K] is the structural stiffness matrix which is 
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built up from the element stiffness matrices 
by a procedure expressible algebraically as 

[K] - [a T ][k][a] 

where [a] is a boolean operator sometimes 
called the connectivity matrix, and [k] 
denotes the element stiffnesses. 


At this point it is appropriate to go through the algebraic 
development of [k] in a step-by-step fashion. 


The displacements within an element vary as defined 
by Equation (4). By differentiating Eq. (4), the required 
tenas are obtained foT substitution in the strain-displace 
ment relationships (2) so that all six strains can be ex- 
pressed in the form 


e l 



e 2 


l b 2 C 2 J 



( 6 ) 


* ' J>12 

'* •* ,ij b 

k 2 ” 1^2 e 2 J 

T - L d 12 e 12 j 



(7) 


where corresponds to the u and v displacements and 

to w. 

Using Eq. (1) for the strain energy expression we can 


U ’ I V 





write 
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where 

tW - f CD m f^I + ‘>I>Lh +b 2-‘-C 1 -' , )Ct»’;} L b 2 J + {bT> L b lJ ) 

* A 

* |-(l-v){bj 2 ) |b 12 J]+D b [{dJ+d|} t d 1+ d 2 i 
- (1-v) ({dj} t d 2 j +{dp udjj )*2(l-v) CdJ 2 ) v d 12 j ])dA 


( 9 ) 

^mb 1 ” f A (D m f{b l* b 2 }lc l +c 2J* (1 ' v)({b I }vC 2 J * {b 2 } lCiJ) 

* A 

* ■y(l" v ) ^12^ l c 12J ^ + ®b^^l +< *2^ L e i +e 2J 


- (l-v)({dJ}Le 2 i*{d^} l e 1 j)+2Cl-vD{dJ 2 ) L e 12 j])dA 


(10) 

I k bmJ * t^ b ] (ID 

and 

[k^fe] * | A CD m [tc 1 +c 2 H c i +c 2J “ ( 1 * v ) l.c 2 j +{c 2 ) (C^j ^ 

♦ 7 (l-v) {c 12 ) l .c 1 2j ] + D^[{e 1 + e 2 ) ie 1 +e 2 J 
- (1-v) ({ej> L e 2 j +{e 2 ) te ) +2(l-v) {e^} L e 12 J ])dA 


where D 

m 


Et 

(1-V 2 ) 


and * 


Et 

12(1 -v 2 ) 


( 12 ) 


Thus the element stiffness matrix [k] is defined. 

The surface integral in Eq. (5) can be expressed in 
terms of nodal displacements by introduction of Eq. (4). 

| pwdA * [Wj | p{N}dA - ^wj {P e > 

A A 
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The element loads {P }, usually referred to as work 

w 

equivalent loads, can be added to the load vector {P} to 
give a total load vector { P^. } : 

{P T } - {P} ♦ [a T ]{P e } 

Equation (5) is rewritten 

n - | L Aj[K]{i} - L*j{ p T } (13) 

The constraint equations are treated next. These equa- 
tions can be expressed as 

[C] {A} =* {0} (14) 


where [C] is an m x n matrix of coefficients which depend 
only on the element geometry (m is the number of constraint 
equations and n the total number of freedoms in the struc- 
ture). {0} is an m x 1 null vector. 


Each component equation in matrix equation (14) is 
derived in the following way. At the midside of an edge 
we impose 


4 


n 


- <J> = 0 

n 


(15) 


where $ n is the rotation of the tangent to the shell which 
is perpendicular to the edge, and r and s are the neighbor- 
ing elements. If 4> n is represented by a vector pointing 
from node i to node j (Fig. 2) , by vector transformation. 


*1 “ -$T s i n 9:; + $^>COS 0.. (16) 

n l ij i ij 

It is more convenient to define as a vector pointing from 
node j to node i. In this way the cyclic ordering of the 
nodal coordinates is utilized for the computation of sin 6^ 
and cos 6^. But because of this the sign in Eq. (15) has 
to be changed, i.e. + *n * ® (15). It so happens that 
this simplifies the algorithm because all terms are added 
into the constraint equation regardless of which element is 
referred to. 

Using the same approach as used in setting up Eqs. (3) 
through (12), by reference to Eq. (2), 
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(18) 

T X T T 

where (£ 2)1 (^ 2 ^ and ^ 82 ^ consist of constant terms ob- 

tained by evaluating the fourth and fifth of Eqs. (2) for the 
midside coordinates. The same procedure is performed to ob- 
tain an equation that corresponds to Eq. (18) for element s. 
Thus all the coefficients for one row of Eq. (14) are estab- 
lished. 

Having established the procedure for deriving the con- 
straint equations, the treatment of these equations is dis- 
cussed next. 

Using the Lagrangian multiplier technique, the poten- 
tial energy expression (13) can be augmented to give 

JI 1 - | l Aj[K](A} - L Aj{P T J + ,A|[C]U} (19) 

Variations of II* with respect to A^ and the Lagrangian 
multiplier can be taken independently. Thus 

|2- - 0 : [K] (A) - (P T ) + [C T ](X} - {0} (20) 

|^ - 0 : [C] { A} - (0) (21) 

Combining Eqs. (20) and (21), 



( 22 ) 
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By treating the constraints as outlined above, the size of 
the problem has been expanded somewhat. Alternatively, the 
number of equations could have been reduced by using Eq. 
(14) to eliminate certain freedoms in Eq. (13). This ap- 
proach is simple but proves to be awkward to incorporate 
in a general purpose program. 

The system of equations represented by (22) are banded 
if the displacements and Lagrangian multipliers are suit- 
ably ordered. The ordering of the freedoms is important 
in one other respect. If a constraint equation appears 
before any of the displacement freedoms to which it is 
coupled, a zero pivot term will result in the Gauss elimina 
tion mode of solution. Pivoting, or careful numbering of 
the degrees-of -freedom circumvents this problem. 


It is worthwhile to point out that the form of Eq. 

(22) is not a handicap in the building of the system stiff- 
ness matrix. The constraint equations can be included in 
the element 'stiffness 1 matrix to give it the following 
structure : « 


k 

c 



( 22 ) 


where [c] is a 3 x 30 matrix, each row corresponding to the 
coefficients given in Eq. (18). Thus for element r, no con- 
straint coefficients from element s would be present in the 
stiffness matrix. This element stiffness can be treated in 
the usual way. 


In some problems in structural analysis, constraint 
equations unrelated to particular elements are imposed. 

These can be treated simply by using the approach outlined 
above if one uses the 'constraint element' device, as pointed 
out in Reference 35. One imagines a nonphysical element 
which connects all the nodes in the constraint equation to 
a 'node' which is associated with the Lagrangian multiplier. 
The corresponding 'stiffness' matrix has the form 



r 

j 
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where c consists of the coefficients in the given constraint 
equation. Thus, general purpose programs that permit the 
direct input of stiffness matrices (e.g. STRUDL) can cope 
indirectly with constraint equations. 


The presence of the centroidal node of a triangular 
element is awkward for the input of data. The degrees of 
freedom at the centroid can be eliminated conveniently by 
means of static condensation. The original element stiff* 
ness matrix (33x33) is partitioned in the following manner. 



k bb 

k bs 


k sb 

k ss 

The subscripts b and s 

refer 



1 


(23) 


refer to the freedoms that are to be 
retained and eliminated respectively. [k^l is (30x30) and 
k ss is (3x3). Solving the two matrix equations of (23) for 

<V 


lk bb ' k bs k ss lfc sb 


HA^} - {F. 


' k bs k ss 


-1 


V 


(24) 


Thus the new modified (30x30) stiffness matrix for the tri- 
angle with no centroid, is the first term of Equation (24). 
The element loadings have also to be modified in the manner 
indicated by the right hand side of Equation (24) . 

For rectangular problems it is advantageous to combine 
triangular elements to form quadrilaterals. Not only does 
this artifice give better computational results, but it 
also yields a more efficient element. For a rectangular 
mesh a saving of almost 50% in the total number of freedoms 
is realized by using the quadrilateral, and there is a cor- 
responding reduction in the bandwidth of the resulting sys- 
tem equations. 

I 
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The method of static condensation is again used to derive 
the quadrilateral element stiffness matrix (40x40) from the 
component triangular element matrices (33x33) . The geometry 
of the quadrilateral is shown in Figure 4. 

VII. NUMERICAL RESULTS 

Numerical results are given in this section for four 
problems, representing a broad range of circumstances en- 
countered in thin shell analysis. The problems are those for 
which alternative numerical solutions, and in some cases clas- 
sical solutions are available and comparisons of the respec- 
tive solutions are presented as a function of parameters which 
define the solution effort. 

The first problem solved is that of a pinched cylinder 
(Figure 5). A classical solution is available for this ease, 
but for a state of deformation that is less general than that 
represented by finite element analysis. The second problem 
is the gravity loaded cylindrical shell of Figure 6. This 
problem is effectively the "standard" basis of comparison for 
finite element shell formulations. Analyses of the hypar shell 
(Figure 7) serve to cover the case of a structure with negative 
Gaussian curvature. The pear-shaped cylinder (Figure 8) is a 
departure from geometric regularity and has been devised 
specifically (Ref. 9) as a test-case for finite element thin 
shell formulations. 

Pinched Cylinder 

This problem, shown in Figure 5a involves diametrically 
opposed point loads P = 100 lb. acting on an isotropic thin 
cylinder of elastic modulus E = 10.5 x 10° lb. /in. , Poisson's 
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ratio v = 0.3125, with length L = 10.35 in. and radius R = 

4.953 in. The problem is widely used for evaluation of new 
shell finite elements because it emphasizes the significance 
of rigid body effects, which many newly -formulated shell ele- 
ments either attempt to approximate or represent exactly. An 
analytical solution for the cylinder with free ends and loaded 
in this manner, but neglecting extensional behavior, is avail- 
able in Ref. 36; the value of the displacement under the point 
load is too small due to the inextensional assumption. There 
is no available classical solution which accounts for extensi- 
bility. The converged finite element solution given by Cantin 
(Ref. 30) is therefore taken to be the correct solution. 

The boundary conditions for finite element analysis con- 
sist of the application of symmetry conditions at the circum- 
ference passing through the point loads and also along the 
longitudinals which divide the shell into 90° segments. (The 
finite element representation therefore encompasses one-eighth 
of the shell) . 

Three different finite element solutions for the displace- 
ment under the load are shown in Figures 5b and 5c as a function 
of the number of solution unknowns. Figure 5b refers to a 
situation wherein each element spans the half-length of cylinder 
grid refinement takes place entirely in the circumferential 
direction. This means that in Fig. 5c there are on the order 
of twice as many degrees-of -freedom for a given level of grid 
refinement in the direction as in the scheme associated with 
Figure 5b. 
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One of the represented finite element solutions is due to 
Cantin (Ref. 30) who employs a rectangular cylindrical element 
based on radial and membrane displacement fields which are 
coupled in such a way that the rigid body motion condition is 
satisfied. Another of the finite element representations shown 
in Figs. 5b and 5c is due to Bogner, et al (Ref. 37). This is 
also the restricted case of a rectangular cylindrical element 
but is based on uncoupled bicubic expansions for the three 
displacement components, u, v and w. It consequently represents 
’'implicit" satisfaction of the rigid body motion condition. 

The third set of results derives from the present formulation. 

Study of Figs. 5b and 5c discloses that results of the 
present formulation lie between those of Cantin and Bogner, et 
al , and are of reasonable accuracy for a modest number of un- 
knowns. Again noting that this problem is very sensitive to 
the rigid body motion condition, it is reassuring that the 
results from the present formulation are satisfactory even 
though the rigid body modes are not explicitly included in the 
displacement functions. 

Comparison of Figs. 5b and 5c discloses that fewer degrees- 
of-freedom are needed for a given level of solution accuracy 
for the "coarse" representation, than with two longitudinal 
elements. This is due to the fact that refinement of grid in 
the longitudinal direction does not improve solution accuracy 
to an extent that counterbalance the large expansion of solution 
unknowns . 

Cylindrical Roof 

This problem concerns a cylindrical shell roof (see Figure 
6) supported by rigid diaphragms at the ends and with free edges , 
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subjected to gravity loading. This particular problem has 
received a great deal of attention from thin shell analysts for 
comparison purposes, since an alternative solution, stemming 
from work described in Ref. 38, is available. 

The cylindrical shell if 50 ft. long, 3.0 in. thick and of 
40 ft. radius. It is composed of isotropic material with 
E = 3 x 10° lb/in ^ and v * 0. Due to the two axes of symmetry 
the analysis is performed for a square grid of elements in one 
quadrant of the shell. 

Figure 6 shows the vertical displacement at the midpoint (A) 
of the free edge as a function of the number of analysis unknowns 
(degrees -of -freedom) . Three finite element solutions are compared. 
One is based on the use of flat plate elements and is taken from 
Ref. 39. The particular elements employed were the constant 
strain triangle for membrane action and the "HCT" triangle (Ref. 

40) for bending. It is seen that the numerical accuracy of these 
results is poor. Carr (Ref. 41) has demonstrated, however, that 
by use of improved flat plate elements it is possible to obtain 
significantly better results than are shown here. Indeed, 
higher-order membrane representation is indicated as a means 
of improving the accuracy of solution. Although flat element 
formulations do not use curved shell strain-displacement equations 
explicitly, they do converge to the deep shell solution. 

The second finite element solution shown in Figure 6 is due to 
Cowper, et al (Ref. 13). This formulation, which represents one 
of the most sophisticated developments using a shallow shell theory, 
was described in Section II. The results are seen to approach 
the analytical shallow shell solution. These same authors discuss, 
in Reference 14, means of correcting their formulation to deal with 
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deep shell behavior; no numerical results are available, how- 
ever, for the present problem. 

The third finite element solution is that given by the 
element formulated in this report. Figure 6 discloses that 
the present results are quite accurate at all levels of grid 
refinement, but that they do not converge monotonically to 
the deep shell solution. It is felt that these results are 
artificially improved for the coarser idealizations because 
of the way the distributed loading has been treated. The load- 
ing is of the gravity type and thus has a sine and cosine 
distribution. The computer program used for the present analy- 
sis accommodated only a linear distribution of pressure between 
the nodal values within each element. Thus the assumed load- 
ing is everywhere slightly smaller than the actual loading, 
but this discrepancy reduces with increasing grid refinement. 
Had the correct variation of pressure been incorporated the 
result for the coarsest grid would have been higher and mono- 
tonic convergence is anticipated. 

Clamped Hypar 

The clamped hypar (Fig* 7) is another problem used exten- 
sively for comparison purposes. The shell has a square plan- 
form and is clamped on all four edges. The loading is a con- 
stant normal pressure p over the complete shell. The ratio 
of the maximum rise of the surface to the horizontal projection 
of the shell edge is 1:5. The ratio of the maximum rise to the 
thickness is 1:25. The material has a Poisson’s ratio of 0.4. 
The shell exhibits double symmetry about the diagonals. Advan- 
tage of this symmetry can be taken by triangular idealizations, 
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whereas the whole shell has to be analyzed by rectangular ele- 
ments. Thus two comparison plots are shown, one for triangular 
elements for a quarter of the shell, and one for all elements 
where the complete structure is analyzed. 

The present element is compared with Ford's curved tri- 
angular element (Ref. 42) in the first plot of Figure 7. Here 
the normal displacement at the shell center is plotted in non- 
dimensional form against the total number of freedoms in the 
idealization. Ford's element converges to the true solution 
at a slightly more rapid rate than the present element. The 
reason for this is the adoption of a more complex quintic dis- 
placement field for the normal component, as compared to the 
cubic approximation adopted herein. 

In the second plot of Figure 7, the results of the present 
study are compared with those of Ford (Ref. 42) and Connor and 
Brebbia (Ref. 43). Again Ford's element has the most rapid 
convergence. The Connor and Brebbia element, which also uses 
a cubic normal displacement field, gives similar accuracy to 
that achieved by the present element, but converges from the 
too flexible side of the exact solution. 

It is noticeable that the results obtained with the present 
element for this problem are inferior to those obtained for the 
cylindrical roof. The reason for this is that bending plays a 
more important role in the shallow hypar. Membrane behavior 
is represented more accurately in the present development. 
Pear-shaped Cylinder 

The problem of a pear-shaped cylinder under axial compres- 
sion shown in Figure 8, has been proposed (Ref. 44) as a test 
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problem for the evaluation of shell analysis computer programs. 
The structure was chosen to simulate the general geometry of 
problems encountered in a certain class of space vehicles. 

This is in contrast to the regular geometry of problems to 
be found in the literature of curved thin shell numerical 
analysis . 

The loading on this structure is defined by a uniform 

axial shortening of 2 x 10 ^ in. The material is aluminum 
7 2 

with E = 10 lb. /in. and v = 0.3. Boundary conditions of 
the simple-support type, detailed in Ref. 44, are invoked. 

Care must be taken in the connection of the flat and curved 
sections of the shell. It can be recalled that the order of 
polynomials used in description of membrane behavior of the 
subject element require derivatives as degrees -of -freedom. 

In a flat element, along the juncture line with curved ele- 
ments, this degree -of -freedom is simply the first derivative 
8 v 

) , where v and y are the membrane displacement and coordi- 
dy 

nate axis perpendicular to the juncture line. In the curved 

segment it is + S, where R is the relevant shell radius 
dy K 

of curvature. If only the continuity of the first derivative 
is enforced the solution will be significantly in error. 
Satisfactory results can be obtained by either using the full 
definition of this degree-of -freedom in the curved segment or 
by neglecting entirely the continuity of this degree -of -freedom. 
(The continuity of the u- and v degrees -of -freedom must of 
course be enforced) . The latter option was adopted in this 
work. This aspect of the analysis of the subject problem is 
discussed in some detail in Ref. 45. 
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In the present study the shell has been idealized by 
quadrilateral elements (4 triangles per quadrilateral) , 

(Fig. 4) with approximately 300 degrees of freedom after 
imposition of boundary conditions. The results of the analysis 
are presented in Figure 9 along with results obtained from a 
finite difference (STAGS-Ref. 46) computer program that 
involved a 4 x 40 gridwork with 756 degrees of freedom. The 
results obtained with the present development compare favorably 
with the alternative solution. No determination was made of 
the solution accuracy as the analysis grid was refined. Ref. 

44 gives further numerical studies of this problem. 

VIII CONCLUDING REMARKS 

The formulation of a triangular thin shell finite element 
has been presented in this report. The three components of 
displacement (u, v, w) are each described in terms of cubic 
polynomial functions. The difference in the order of bending 
versus stretching derivatives in the strain energy expression 
implies lower-order membrane (u , v) fields than the bending 
(w) field. The use of equal -order fields in the present work 
has the effect of approximating more closely the rigid body 
motion condition. The interelement displacement continuity 
condition is not met by the basic element displacement fields, 
but this is also approximated more closely by writing the 
continuity requirements as conditions of constraint. These 
constraint conditions are imposed by use of the Lagrange 
multiplier technique. 
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The present formulation was demonstrated to be reliable 
in a very wide range of linear analysis circumstances. Its 
development in terms of deep, rather than shallow, shell 
theory was verified to be of importance. Also, in the test 
case that emphasizes the rigid body motion condition, the 
approximate satisfaction of this condition was demonstrated 
to lead to results of acceptable accuracy. 

The standing of this formulation, when compared with 
alternatives, must take into account the relative simplicity 
of the chosen displacement fields. The question of simplicity 
is also important in the intended application to geometrically 
nonlinear analysis, where substantially more complex expres- 
sions for the stiffness coefficients must be dealt with. The 
extension to nonlinear analysis is taken up in a companion 
report (Ref. 6) . 
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APPENDIX A 

ASSESSMENT OF EFFECT OF VIOLATING CONDITIONS 
ON RIGID BODY MOTION 


Herein an attempt is made to estimate the strain 
energy content of a curved beam element when subjected to 
a rigid body translation. It is assumed that the displace- 
ment fields of that element do not explicitly describe such 
a mode of displacement. 

In Figure 3 a curved beam is shown being displaced 
rigidly by a distance r. The axial and normal displace- 
ments u, and w are 

e 3 e 5 

u * -rsin 0 * -r ( 0 - jy* + ••••) 

2 4 (A-l 

w * rcos 6 *r(l“ 2 T + 3 “j- 


If the mode of representation of both the radial (w) 
and circumferential (ii) displacements is cubic in the 
angular coordinate 0 , we have 


u * a„ 
o 

+ a 1 0 

+ a 2 e 2 

+ a 3 

e 3 

* 

K 

O* 

o 

+ b^0 

+ b 2 e 2 

+ b 3 

e 3 


(A- 2) 


and by comparison with (A-l) we can write the rigid-body-mode 
form of these components as 


u * -r 0 



w • r 



(A-3) 


Substituting these displacement functions into the relevant 
strain energy expression, we find 
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" | 0 ^ (1 ‘ V * I fl * T -51 dx 
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♦ ^[ e '-h + h v- 4)J 2 dx 
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40 

where the subscript 3 signifies the cubic function for u. 

For a linear function for u, can be calculated as 
follows: 

U 1 - ^7 \ 0 C TT + R " IR 3 Rd0 + J G C ^f + ^T } Rd0 


U 


1 


2*5 

EAr^e* 

TOE 


Now, a measure of the significance of employing a 
cubic function for the circumferential displacement is 

the ratio between and (each of which would be zero if 
the condition on zeTO strain under rigid body motion were met) . 



For a beam of rectangular section, depth t, the ratio be* 
comes 




2 



XV 


u 


3 


« U 


1 


Therefore, the strain energy corresponding to a rigid body 
translation of the refined arch element is negligible com* 
pared to that of the basic element, which itself may well 
be small. 


Figure 1. GEOMETRY OF SHELL ELEMENT 


Figure 2 



TRIANGULAR COORDINATES AND COMPATIBILITY OF 
ROTATION BETWEEN NEIGHBORING ELEMENTS 
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Figure 3. RIGID BODY MOTION OF CURVED BEAM 



O Retained nodal point 
X Eliminated nodal point 


Figure 4. GEOMETRY OF QUADRILATERAL ELEMENT 


Norma! Displacement at Pain! Load (in.) 





Analysis Octant 



FIGURE 5. COMPARISON OF CURVED SMELL ELEMENTS FOR CYLINDER LOADED BY TWO DIAMETRICALLY 
OPPOSITE EXTERNAL CONCENTRATED FORCES AT MIDLENGTH. 
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E = 3.0 x I0 6 lb/in 2 
V =0 

t * 3.0 in 
Load = 90 lb /ft 2 



No. of Unknowns 

FIGURE 6. ELEMENT SOLUTIONS FOR CYLINDRICAL 
SHELL ROOF UNDER GRAVITY LOADING. 
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Z,W 



Pressure Loading * p 
x iq 5 Thickness * t 



No. of Unknowns 


FIGURE 7 COMPARISON OF VARIOUS CURVED ELEMENTS 
IN ANALYSIS OF A UNIFORMLY LOADED CLAMPED 
HYPERBOLIC PARABALOID. 




Circumferential Angle (deg) 0 

Figure 9. COMPARISON OF NORMAL DISPLACEMENT SOLUTIONS ALONG THE LOADED EDGE 
FOR PEAR-SHAPED CYLINDER 





